Phylogenetic analyses reveal bat communities in Northwestern Mexico harbor a high diversity of novel cryptic ectoparasite species

Abstract Parasites are integral parts of ecosystem function and important drivers of evolutionary processes. Characterizing ectoparasite diversity is fundamental to studies of host–parasite interactions, evolution, and conservation, and also for understanding emerging disease threats for some vector borne pathogens. With more than 1400 species, bats represent the second most speciose mammalian clade, but their ectoparasite fauna are poorly known for most species. We sequenced mitochondrial Cytochrome Oxidase C subunit I and nuclear 18S ribosomal gene fragments, and used Bayesian phylogenetic analyses to characterize ectoparasite taxon identity and diversity for 17 species of parasitized bats sampled along the Baja California peninsula and in Northwestern Mexico. The sequence data revealed multiple novel lineages of bat bugs (Cimicidae), flies (Nycteribiidae and Streblidae), and ticks (Argasidae). Within families, the new linages showed more than 10% sequence divergence, which is consistent with separation at least at the species level. Both families of bat flies showed host specificity, particularly on Myotis species. We also identified new records for the Baja peninsula of one tick (Carios kelleyi), and of five Streblid bat fly species. One Nycteribiid bat fly haplotype from Pallid bat (Antrozous pallidus) hosts was found throughout the peninsula, suggesting potential long distance co‐dispersal with hosts. Different bat bug and tick communities were found in the north and south of the peninsula. This study is the first systematic survey of bat ectoparasites in the Baja California peninsula, revealing novel lineages that are highly genetically differentiated from other parts of North America. For some ectoparasite species, haplotype distributions may reflect patterns of bat migration. This work is a first step in characterizing ectoparasite diversity over the Baja California peninsula, and understanding how ecological and evolutionary interactions shape bat ectoparasite communities among host species in different parts of their ranges.


| INTRODUC TI ON
Bat flies comprise two main families, Nycteribiidae and Streblidae, which have a common origin from a single lineage coevolving with bats (Dittmar et al., 2006;Wenzel & Tipton, 1966). As of 2006, there were approximately 520 species described (Dittmar et al., 2006). With new records added regularly (Saldaña-Vázquez et al., 2019;Szentiványi et al., 2019), taxonomy at the species level is often ambiguous, and an updated global review and synthesis is yet to be produced.
Ticks are grouped in three families: Argasidae (soft ticks), Ixodidae (hard ticks), and Nuttalliellidae. The latter family comprises a single known species, Nuttalliella namaqua (Burger et al., 2014;Mans et al., 2012), which has been suggested as the basal tick lineage (Mans et al., 2012). There are approximately 894 known species of ticks worldwide, with 32 species of Argasidae and 68 species of Ixodidae in Mexico Pérez et al., 2014).
Bat tick phylogenetic studies have mainly focused on old world species (Hornok et al., 2017), with a limited number for the Americas and other parts of the world (Black et al., 1997;Burger et al., 2014), and suggest bat-associated ticks commonly exhibit host-specificity (Sándor et al., 2019). Classification of taxonomic relationships among soft ticks are controversial (Burger et al., 2014), and more studies are needed to accurately describe their taxonomic status (Burger et al., 2014;Estrada-Pena et al., 2010).
Baja California (Figure 2) is an isolated peninsula, 1300 km long, in northwestern Mexico. Its complex geological history means it has a mosaic of habitats, with high levels of biodiversity and endemism K E Y W O R D S 18S, Argasidae, Chiroptera, Cimicidae, COI, DNA barcoding, ectoparasites, Nycteribiidae, phylogeny, Streblidae

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology, Ecological genetics, Genetics, Parasitology, Phylogenetics, Taxonomy F I G U R E 1 Specimen of canyon bat (Parastrellus hesperus) infested by bat bugs, sampled for this study in northwestern Mexico. Image credits: The authors, University of Leeds.   (González-Abraham et al., 2010;Nájera-Cortazar et al., 2015;Riddle et al., 2000), with approximately 25 species of bats recorded Medellín et al., 2007). While there are several studies describing Baja bat diversity, distributions, and ecology Frick et al., 2008;Nájera-Cortazar et al., 2015), there are no studies describing the bat ectoparasite fauna.
Final extension was performed at 72°C for 10 min.
PCR products were visualized by gel electrophoresis using 1% agarose with GelRed® (Biotium Ltd) staining, and then sent to Genewiz Inc. (Azenta Life Sciences), for PCR product purification and Sanger sequencing, with each amplicon sequenced in both forward and reverse directions.

| Quality control and reference sequences
Sequence quality for both the forward and reverse strands of each amplicon was evaluated in BioEdit 7.2.5 (Hall, 2005). Trimmed forward and reverse sequences were combined to generate a consensus sequence for each amplicon, and then analyzed in BLAST (The National Library of Medicine, 2018) to generate initial taxon identities and identify reference sequences. Additionally, COI sequences were also analyzed in the BOLD system platform (Ratnasingham & Hebert, 2007).
Reference sequences for phylogenetic analyses were compiled from previous studies on each ectoparasite group (Burger et al., 2014;Dittmar et al., 2006;Mans et al., 2012;Tortosa et al., 2013), and by performing a systematic search using AnnotationBustR 1.2 package (Borstein, 2018) in RStudio 1.1.456 (RStudio Team, 2021), searching for the closest genus for the sequences generated in this study and for those obtained by BLAST. Alignments combining reference sequences and those from this study were generated using CLUSTAL W (Thompson et al., 1994), implemented in BioEdit, and were reviewed by eye, with manual correction of potential errors as required.

| Sequence summary statistics and phylogenetic analyses
Genetic distance estimates among haplotypes and best fit sequence evolution models for the COI and 18S datasets were evaluated using MEGA 10.1.7 (Sudhir et al., 2018). The Barcode Index Number system was followed to delimit a lineage, where intraspecific variation at COI is generally considered as groups of sequences with less than 2% divergence, and exhibiting more than 4% divergence from neighboring lineages (Hebert et al., 2003;Ratnasingham & Hebert, 2013;Salinas-Ramos et al., 2018).
Bayesian phylogenetic analyses were performed using the best fit evolution model identified for each group and each marker implemented with BEAST 1.10.4 . A MCMC chain length of 10,000,000 was used and priors specific to each parasite group and sequence evolution model were selected using the program BEAUti 1.10.4 . Phylogenetic analyses were performed separately for each ectoparasite group and each marker. Nycteribiidae and Streblidae families of bat flies were analyzed separately in order to improve sequence alignment quality and phylogenetic resolution. For each family-marker combination, two separate runs using the same Bayesian settings file generated by BEAUti were run in BEAST, with a burn in of the 10% of the total number of iterations. After, stationarity of BEAST results were assessed in Tracer 1.7.1 . Both files for each ectoparasite family-marker combination were combined using LogCombiner 1.10.4 , generating a single .log file and a single .tree file. A majority-rule consensus tree was inferred using TreeAnnotator 1.10.4  with the combined .tree file, and using a posterior probability limit of 0.6, median nodes heights were summarized. Phylogenetic trees were annotated, edited, and visualized using iTOL 5.6.2 (Letunic & Bork, 2016).
Following phylogenetic analysis, haplotypes for each family were grouped by lineage, and haplotype diversity and genetic summary statistics were calculated in DNAsp 5.10.01 (Librado & Rozas, 2009   Flies were the predominant ectoparasite type, found in 13 of 17 sites, followed by ticks and bugs, present in eight sites each. Initial field morphological evaluations suggested most bug specimens were related to Cimex pilosellus (Usinger, 1966). Most bat flies could be identified to the genus level, which were later corroborated with molecular data. On the basis of morphology, all bat ticks were identified as family Argasidae (soft ticks), with at least six different morphotypes present. Most of these were tentatively attributed to the genus Ornithodoros. there were no species level matches for neither COI nor 18S markers.

| BLAST and BOLD sequence matches
In the case of streblid bat flies, 14 sequences matched (~99%-100% in GenBank-BOLD databases) with five known species for COI (see Table 2); but there were no matches for 18S. These marker differences are likely attributable to database coverage, since where reference sequences were available for COI and 18S from the same specimen/ study, BLAST results for our sequences were consistent for both markers. For brevity and to provide comparability with larger numbers of reference sequences further reporting of diversity, divergence statistics, and taxonomic identity will be focused on COI results.

| Genetic diversity
Genetic diversity statistics for each COI lineage in each ectoparasite family are summarized in Table 2. Excluding four lineages represented by single individuals (e.g., Cimex 1), and three lineages including only individuals with the same haplotype (e.g., Basilia 1), nucleotide diversity (N d ) ranged from 0.001 to 0.006, with the highest value 0.071, presented by the Tick 4 lineage, with three haplotypes (H) in three sequences ( Table 2). In general, the number of haplotypes were close to or the same as the number of sequences tested per lineage.

| Phylogenetic assessment of Baja peninsula bat bug sequences
The best fit evolution model for the COI gene set was GTR + G + I, and K2 + G + I for 18S. For each marker, 30 sequences were generated, forming four novel lineages with respect to reference sequences ( Figure 4). Genetic divergence between the four peninsular lineages ranged from 9.9% to 17.1% (Appendix 1), and between 7.2% and 20.9% against reference sequences, with C. latipennis, C.
antennatus, and C. adjuntus presenting the lowest divergence values.
Phylogenetically, all the lineages sit within the Pilosellus complex, of North American members of the Cimex genus.
For COI, Cimex 1 is represented by a single specimen (EBCO155) obtained from a Myotis yumanensis host at Mosqueda, and forms a sister lineage to Cimex 2, represented by two haplotypes from two bugs parasitizing M. californicus hosts, which were also sampled from northern sites. The 18S analysis also placed EB18S155 in a separate clade, but Cimex 2 is paraphyletic. Both markers placed Cimex latipennis, which is distributed from Canada to the north western USA (Usinger, 1966), as the closest named molecular reference species to these 2 lineages. Cimex 3 is represented by a monophyletic clade for COI, consisting of two haplotypes, sampled from specimens col-

| Phylogenetic assessment of bat fly sequences
The best fit sequence evolution models were GTR + G for the COI gene set, and T92 + G + I for the 18S data. In total, 77 bat flies were  (1966)
Supporting this, it is noted that each lineage had different host species, Sturnira parvidens and Artibeus hirsutus, respectively ( Table 2).
Most streblid lineages were parasitizing fruit-nectar feeding bats Posterior probability support is indicated by the size of black circles at tree nodes. Where available, information for location and host is written next to each reference sequence label. To improve clarity of the tree, collapsed clades of reference sequences are shown as gray triangles. GenBank accession numbers for reference sequences are given in the sequence labels. and Primavera, Figure 2, Table 1). Outside of the closest sequence matches, none of the other species reported as parasitizing hosts in Baja (Table 2) appear to have reference sequences in GenBank.

Genetic divergence within groups from the lineages in this study
showed levels between 0.0% and 1.1%, with the exception of A.

| Phylogenetic assessment of Baja peninsula bat tick sequences
The best fit sequence evolution models were GTR + G + I for the COI gene set, and K2 + G + I for the 18S gene set. There were 45 ticks sequenced from the Argasidae family (soft ticks), with 39 sequences generated for COI, and 44 for 18S, where six specimens failed to amplify for COI (specimens ETCO157, ETCO306, ETCO338, ETCO480b, ETCO480c, and ETCO485), and one specimen failed to amplify for 18S (specimen ET18S457). One new Baja species record and seven potential novel lineages were obtained (Figure 7).
The only match with GenBank and BOLD COI records (95.95% and 97.07%, respectively) was the soft tick Carios kelleyi, for 11 specimens (parasitizing A. pallidus hosts) with intra-clade divergence of 1.94%-2.47% (Figure 7, Appendix 4) yumatensis showed divergence of around 9.9% with respect to Tick F I G U R E 5 Bayesian phylogenetic trees for bat flies of the family Nycteribiidae obtained using mitochondrial COI (left) and the ribosomal 18S (right) sequences. Posterior probability support is indicated by the size of black circles at tree nodes. Where available, information on location and host is written next to each reference sequence label. To improve clarity of the tree, collapsed clades of reference sequences are shown as gray triangles. GenBank accession numbers for reference sequences are given in the sequence labels.
6 lineage. Ticks of lineage Tick 6 were recovered from M. peninsularis from which there are no previous records of bat ticks to our knowledge.
All sequences from this study generated the same lineages for both COI and 18S genes (Figure 7). Lineages C. kelleyi and Tick 1 to 5 formed a monophyletic group with respect to the reference sequences with posterior support greater than 0.6 for both markers, while Antricola 1 and Tick 6 lineages were positioned in separate clades. The topology of sister clade relationships varied slightly between markers, particularly around deeper nodes which had posterior probability support less than 0.85. This suggests more data is required to resolve deeper taxonomic relationships among species.
Antricola 1 and Tick 6 were separated around shallow deep nodes in both analyses, where the absence or presence of reference sequences influenced their topological proximity (Figure 7). None of the reference sequences from Antricola, Carios, and Ornithodoros genera, form monophyletic groupings with respect to genus nomenclature.
Tick lineages 2 and 6 (gray and orange, Figure 7) were only ob- marginatus reference was available. The two host species for ticks of this lineage were, sampled mostly at mid-peninsula (Figures 2 and   3). Tick 5 lineage had the most diverse host and spatial distribution, being found on four bat species, and at sites from the north to the south of the peninsula (yellow clade, Figure 7).  (Hebert et al., 2003), indicating Cimex 1, 3, and 4, could be classed as novel species under this criterion. Genetic divergence of Cimex 2 compared with C. latippenis was 3.2%, but C. latippenis has not previously been recorded parasitizing M. yumanensis (Braun et al., 2015), which might also suggest Cimex 2 as a potential new cryptic species.

| DISCUSS ION
Cimicid bugs have low inherent dispersal capacity, generally feeding for a few days, before dropping from the bat host to digest the blood in the roost, where they can survive without feeding for approximately 1.5 years (Ossa et al., 2019). The population structure of bat bugs is mainly influenced by bat movements (Ossa et al., 2019;Talbot et al., 2016Talbot et al., , 2017Usinger, 1966). While the current data in-

| Bat flies
Ten novel genetic lineages and six new records of bat flies for the study region were found. Nycteribiid flies were more abundant in the northern temperate sites, while streblids were more abundant in the southern and subtropical sites, supporting trends noted by Dittmar et al. (2006). To our knowledge, molecular records of Basilia In previous studies it was found that host-specificity can vary according to the diversity and geographic distribution of hosts (de Vasconcelos et al., 2016;Graciolli et al., 2007;Saldaña-Vázquez et al., 2019). We found that most of our Nycteribiid lineages exhibited host specificity, despite having hosts with overlapping ranges and which may share roost sites (e.g., Myotis species in northern Baja). The Nycteribiid 2 lineage was found on multiple Myotis species, which could potentially indicate specificity at a genus level. Streblid winged flies have previously been described as mostly non host-specific (Dittmar et al., 2006), but more recent studies suggest most species to be host specific ( M. vivesi is endemic to the Gulf of Cortes, and restricted to coastal habitats because of its piscivorous diet (Blood & Clark, 1998;Herrera-Montalvo et al., 2017). A. pallidus is sympatric with M.
vivesi on the Gulf coast but has a wider distribution across western North America. It is primarily an insectivorous feeder, but also includes scorpions and nectar in its diet (Frick et al., 2009). Basilia fly records previously reported for M. vivesi, but without reference sequences, include B. plaumanni, B. pynzonix, and B. producto (Graciolli et al., 2007), while flies parasitizing A. pallidus have been described as B. antrozoi ( Table 2). The threshold level of divergence between the Basilia 2a and Basilia 2b lineages may indicate recent divergence from a common ancestor, and potential incipient speciation driven by association with sympatric but ecologically differentiated hosts.

| Bat ticks
A new record of Carios kelleyi and seven novel tick lineages belonging to the Argasidae family were found. For both COI and 18S tick lineages 1-5 and Carios kelleyi formed a clade with more than 0.6 posterior support, suggesting they form a Baja or western North America endemic lineage of bat ticks derived from a common ancestor.
The genera Antricola, Carios, Nothoaspis, and Ornithodoros are associated with bats and their roosting sites in Mexico (Sánchez-Montes et al., 2016). Compared with Ixodidae, the family Argasidae has few molecular studies and reference sequences (Porter & Hajibabaei, 2018). Classifications at genus level are controversial, with many genera being paraphyletic in existing phylogenies, and debates over synonymous use of genus names such as Carios and Ornithodros (Burger et al., 2014). Such

F I G U R E 6
Bayesian phylogenetic trees for bat flies of the family Streblidae obtained using mitochondrial COI (left) and the ribosomal 18S (right) sequences. Posterior probability support is indicated by the size of black circles at tree nodes. Where available, information on location and host is written next to each reference sequence label. To improve clarity of the tree, collapsed clades of reference sequences are shown as gray triangles. Genbank accession numbers for reference sequences are given in the sequence labels.
ambiguities are also reflected in our phylogenies. Furthermore, while many key internal nodes had posterior support greater than 0.6, differences in reference sequence availability made it challenging to interpret the consistency of inter-clade placement between makers.
Argasid ticks have previously been reported for the bat species in this study, primarily Ornithodoros species (Table 2), but with the exception of Carios (Ornithodros) kelleyi, none are close molecular matches for our sequences. In ticks, for COI, the threshold for between genus divergence is considered to be above 10% (Hebert et al., 2003). The observed COI divergence among lineages of this study (6.1% up to 19.3%), and to reference sequences (9.9%-22.8%), suggests that the Baja lineages could represent novel species and potentially novel genera.
We found apparent host-specific and generalist lineages for the ticks reported here. Lineages that appeared host-specific were restricted to single sites, while generalists were found across multiple sampling locations. For example, Tick 3 parasitizing M. yumanensis F I G U R E 7 Bayesian phylogenetic trees for bat ticks of the family Argasidae obtained using mitochondrial COI (left) and the ribosomal 18S (right) sequences. Posterior probability support is indicated by the size of black circles at tree nodes. Where available, information on location and host is written next to each reference sequence label. To improve clarity of the tree, collapsed clades of reference sequences are shown as gray triangles. GenBank accession numbers for reference sequences are given in the sequence labels.
was found only in San Basilio, and Tick 6 was found only on M. peninsularis sampled in Teste. Carios kelleyi was found to only parasitise A. pallidus in this study, and was found at three sites in the mid-and north peninsula. However, C. kelleyi is known to parasitise multiple species across North America and Cuba (Gill et al., 2004), and the reference COI sequence used here derives from an Eptesicus fuscus bat sampled in New Jersey, eastern USA (GenBank accession code: MT780277). Argasid ticks show a continuum of hosts-specificity (Cumming, 1998;Esser et al., 2016), but tick stage-cycle must be considered, with immature ticks being more generalist than adult conspecifics (Esser et al., 2016;Nava & Guglielmone, 2013). In this study, life stage was not assessed while collecting ticks; therefore, it is possible that there are gaps in host range regarding unidentified larvae that were not sequenced.

| Limits on ectoparasite sampling and identification
The present study identified 21 novel genetic lineages, plus 7 new ectoparasite species records, from 138 bats of 17 species, sampled across 18 sites and 2 years. This suggests a diverse ectoparasite fauna in this previously unsurveyed region of Mexico, but it is also likely to be an underestimate of the true diversity, due to constraints around sampling effort. Bat sampling was limited to May-August and conducted at relatively accessible locations with water sources to facilitate bat capture. While our sampling sites were chosen to be representative of habitat types across Baja, increasing the spatial and temporal scope of sampling would be likely to increase the number of ectoparasite species discovered. Assessment of ectoparasite fauna found in roosting sites against those found feeding directly from their hosts and expanding seasonal coverage will be important for future work.
Although previous studies report limited data on bat ectoparasites from North Western Mexico and South-Western USA (Bradshaw & Ross, 1961;Braun et al., 2015;Pérez et al., 2014;Usinger, 1966), they do not integrate morphological and molecular information. For many species, no reference sequence is available from voucher specimens, or there are errors in species identifications and incorrect annotation of reference sequences. Therefore, for all the ectoparasite groups in this study, further work is needed to unify molecular and morphological characterization, to fully confirm which lineages represent previously undescribed species, and which are species with morphological descriptions but no previous molecular record. This is particularly important for bat tick lineages where input from expert morphologists is needed to account for different life stages.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

A PPE N D I X 3
Estimates of percentage divergence (COI) over sequence-pairs between groups (diagonal left matrix), and within groups (rightmost column) of bat flies (family Streblidae), showing the lineages obtained in this study (row label shaded) versus closely related reference sequences from GenBank. Names in the first column have the genus abbreviated. n/c -not calculated as only a single sequence available.

A PPEN D I X 4
Estimates of evolutionary divergence (COI) in percentage over sequence-pairs between groups (diagonal left matrix), and within groups (rightmost column) of bat ticks (family Argasidae), showing the lineages obtained in this study (row label shaded) versus closely related reference sequences from GenBank. n/c -not calculated as only a single sequence available.